rm(list=ls())

pacman::p_load(lubridate, tidyverse, haven, gtrendsR, rdrobust, rddtools,
               broom, ggpubr, gridExtra, dmetar, meta, zipcodeR)

setwd('put_your_wd_here')

runRdit <- function(out, name){
  sd_est <- sd(out$n,na.rm=T)
  house_rdd <- rdd_data(y=out$n, x=out$running_day, cutpoint=0)
  bw_ik <- rdd_bw_ik(house_rdd)
  reg_nonpara <- rdd_reg_np(rdd_object=house_rdd, bw=bw_ik)
  fullest <- data.frame(Coef=reg_nonpara$coefficients/sd_est,
                        lower=(reg_nonpara$coefficients - 1.96 * reg_nonpara$coefMat[2])/sd_est,
                        upper=(reg_nonpara$coefficients + 1.96 * reg_nonpara$coefMat[2])/sd_est,
                        shooting=name)
  return(fullest)
}

# shootings
allshootings <- read_csv('datasets/list_of_shootings.csv',col_select = 1:4)
allshootings$date <- mdy(allshootings$date)

# geocode
# df <- read_csv('change_petition_walmart.csv')
# df$date <- lubridate::mdy(df$date)
# df <- df %>% filter(country == 'US')
# 
# df$zip <- normalize_zip(df$zip)
# df$zip <- str_sub(df$zip,1,5)
# df$zip <- str_pad(df$zip,5,'left',0)
# 
# for(i in 1:nrow(df)){
#   df$state2[i] <- tryCatch({reverse_zipcode(df$zip[i])$state}, error=function(e){return(NA)})
#   print(i)
# }
# 
# missing_state_index <- which(is.na(df$state2))
# df$state2[missing_state_index] <- df$state[missing_state_index]
# 
# write_csv(df,'datasets/change_petition_walmart_stategeocode.csv')

df <- read_csv('datasets/change_petition_walmart_stategeocode.csv')
df <- df %>%
  group_by(date, state2) %>%
  count() 
dates <- data.frame(date=seq(min(df$date), max(df$date),by='day'))
formerge <- expand.grid(date=dates$date, state2=unique(df$state2))
df_new <- left_join(formerge, df, by=c('date','state2'))
df_new$n[is.na(df_new$n)] <- 0

mindate <- min(df$date) + months(2)
maxdate <- max(df$date) - months(2)
dates <- allshootings %>%
  filter(date >= mindate & date <= maxdate & include==1)

dat_instate <- vector('list', nrow(dates))
dat_outstate <- vector('list', nrow(dates))
states <- vector('list',nrow(dates))
states[[1]] <- c('GA','CO')
states[[2]] <- c('IN')
states[[3]] <- c('CA')

# subset in and out of state datasets
for(i in 1:nrow(dates)){
  dat_instate[[i]] <- df %>%
    filter(date > dates$date[i] - months(2) & date < dates$date[i] + months(2) &
             state2 %in% states[[i]]) %>%
    group_by(date) %>%
    summarise(n=sum(n, na.rm=T)) %>%
    mutate(running_day = as.numeric(date - dates$date[i]))
  dat_outstate[[i]] <- df %>%
    filter(date > dates$date[i] - months(2) & date < dates$date[i] + months(2) &
             !state2 %in% states[[i]]) %>%
    group_by(date) %>%
    summarise(n=sum(n, na.rm=T)) %>%
    mutate(running_day = as.numeric(date - dates$date[i]))

  print(i)
}

model_results_instate <- vector('list', nrow(dates))
model_results_outstate <- vector('list', nrow(dates))

model_results_instate[[1]] <- runRdit(dat_instate[[1]], dates$shooting[[1]])
model_results_outstate[[1]] <- runRdit(dat_outstate[[1]], dates$shooting[[1]])
# model_results_instate[[2]] <- runRdit(dat_instate[[2]], dates$shooting[[2]])
model_results_outstate[[2]] <- runRdit(dat_outstate[[2]], dates$shooting[[2]])
model_results_instate[[3]] <- runRdit(dat_instate[[3]], dates$shooting[[3]])
model_results_outstate[[3]] <- runRdit(dat_outstate[[3]], dates$shooting[[3]])
  
model_results_walmart <- bind_rows(
  bind_rows(model_results_instate) %>% mutate(var = 'In State'),
  bind_rows(model_results_outstate) %>% mutate(var = 'Out State')
  ) %>% 
  mutate(outcome = 'Gun Control Petition', petition = 'Walmart Guns') %>% 
  remove_rownames()

# emma taylor petition data

df <- readxl::read_excel('datasets/petition_signatures_emma taylor.xlsx')
df <- janitor::clean_names(df)
df$date <- ymd(df$signed_on)
df <- df %>%
  select(state, date) %>%
  group_by(state, date) %>%
  count()
date_range <- range(df$date, na.rm=T)
formerge <- expand.grid(date=seq(date_range[1],date_range[2], by='day'), state=unique(df$state))
df <- left_join(formerge, df, by=c('state','date'))
df$n[is.na(df$n)] <- 0

mindate <- min(df$date) + months(2)
maxdate <- max(df$date) - months(2)
dates <- allshootings %>%
  filter(date >= mindate & date <= maxdate & include==1)

dat_instate <- vector('list', nrow(dates))
dat_outstate <- vector('list', nrow(dates))
states <- vector('list',nrow(dates))
states[[1]] <- c('FL')
states[[2]] <- c('NM')
states[[3]] <- c('PA')
states[[4]] <- c('CA')
states[[5]] <- c('IL')
states[[6]] <- c('VA')
states[[7]] <- c('CA','TX','OH')
states[[8]] <- c('GA','CO')
states[[9]] <- c('IN')
states[[10]] <- c('CA')

for(i in 1:nrow(dates)){
  dat_instate[[i]] <- df %>%
    filter(date > dates$date[i] - months(2) & date < dates$date[i] + months(2) &
             state %in% states[[i]]) %>%
    group_by(date) %>%
    summarise(n=sum(n, na.rm=T)) %>%
    mutate(running_day = as.numeric(date - dates$date[i]))
  dat_outstate[[i]] <- df %>%
    filter(date > dates$date[i] - months(2) & date < dates$date[i] + months(2) &
             !state %in% states[[i]]) %>%
    group_by(date) %>%
    summarise(n=sum(n, na.rm=T)) %>%
    mutate(running_day = as.numeric(date - dates$date[i]))
  
  print(i)
}

model_results_instate <- vector('list', nrow(dates))
model_results_outstate <- vector('list', nrow(dates))

#model_results_instate[[1]] <- runRdit(dat_instate[[1]], dates$shooting[[1]])
model_results_outstate[[1]] <- runRdit(dat_outstate[[1]], dates$shooting[[1]])
#model_results_instate[[2]] <- runRdit(dat_instate[[2]], dates$shooting[[2]])
model_results_outstate[[2]] <- runRdit(dat_outstate[[2]], dates$shooting[[2]])
#model_results_instate[[3]] <- runRdit(dat_instate[[3]], dates$shooting[[3]])
#model_results_outstate[[3]] <- runRdit(dat_outstate[[3]], dates$shooting[[3]])
#model_results_instate[[4]] <- runRdit(dat_instate[[4]], dates$shooting[[4]])
#model_results_outstate[[4]] <- runRdit(dat_outstate[[4]], dates$shooting[[4]])
model_results_instate[[5]] <- runRdit(dat_instate[[5]], dates$shooting[[5]])
model_results_outstate[[5]] <- runRdit(dat_outstate[[5]], dates$shooting[[5]])
#model_results_instate[[6]] <- runRdit(dat_instate[[6]], dates$shooting[[6]])
model_results_outstate[[6]] <- runRdit(dat_outstate[[6]], dates$shooting[[6]])
model_results_instate[[7]] <- runRdit(dat_instate[[7]], dates$shooting[[7]])
model_results_outstate[[7]] <- runRdit(dat_outstate[[7]], dates$shooting[[7]])
model_results_instate[[8]] <- runRdit(dat_instate[[8]], dates$shooting[[8]])
model_results_outstate[[8]] <- runRdit(dat_outstate[[8]], dates$shooting[[8]])
#model_results_instate[[9]] <- runRdit(dat_instate[[9]], dates$shooting[[9]])
model_results_outstate[[9]] <- runRdit(dat_outstate[[9]], dates$shooting[[9]])
model_results_instate[[10]] <- runRdit(dat_instate[[10]], dates$shooting[[10]])
model_results_outstate[[10]] <- runRdit(dat_outstate[[10]], dates$shooting[[10]])

model_results_doctor <- bind_rows(
  bind_rows(model_results_instate) %>% mutate(var = 'In State'),
  bind_rows(model_results_outstate) %>% mutate(var = 'Out State')
) %>% 
  mutate(outcome = 'Gun Control Petition', petition = 'Doctor Guns') %>% 
  remove_rownames()

# meta
calcMeta <- function(data, name, outcome, petition){
  m.gen <- meta::metagen(TE = Coef,
                   seTE = se,
                   studlab = shooting,
                   data = data,
                   sm = "SMD",
                   fixed = FALSE,
                   random = TRUE)
  results <- tibble(shooting = name, 
                    Coef = m.gen$TE.random,
                    lower = m.gen$lower.random,
                    upper = m.gen$upper.random,
                    outcome = outcome,
                    petition=petition)
  return(results)
}

# combine
metain <- bind_rows(model_results_doctor, model_results_walmart) %>%
  mutate(se = (upper - Coef)/1.96) %>%
  filter(var == 'In State') %>%
  calcMeta(., 'Meta Estimate (In State)','Gun Control Petition','Meta Analysis')

metaout <- bind_rows(model_results_doctor, model_results_walmart) %>%
  mutate(se = (upper - Coef)/1.96) %>%
  filter(var == 'Out State') %>%
  calcMeta(., 'Meta Estimate (Out State)','Gun Control Petition','Meta Analysis')

# for full image 
datforplot <- expand.grid(
  shooting = unique(allshootings$shooting[allshootings$include==1]),
  var = c('In State','Out State'),
  petition = c('Doctor Guns','Walmart Guns')
) %>%
left_join(bind_rows(model_results_doctor, model_results_walmart), by=c(
  'shooting','var','petition'
)) 
datforplot$shooting <- factor(datforplot$shooting, levels=rev(allshootings$shooting))
datforplot$var[datforplot$var == 'Out State'] <- 'Out of State'
datforplot$petition[datforplot$petition == 'Doctor Guns'] <- 'Physicians for\nGun Control'

# petition_plot <- datforplot %>%
#   ggplot(aes(x=shooting, y=Coef, ymin=lower, ymax=upper, shape=var)) +
#   geom_hline(yintercept=0, color='red',linetype=2) +
#   geom_pointrange(position=position_dodge(width=0.5), fill='white') +
#   coord_flip() + 
#   facet_wrap(~petition) +
#   theme_bw() +
#   theme(legend.position='bottom') +
#   scale_shape_manual(values=c(21,23), name='') +
#   labs(x='', y='')
  
meta_plot <- bind_rows(metain,metaout) %>%
  mutate(shooting = case_when(
    shooting == 'Meta Estimate (In State)' ~ 'In State',
    shooting == 'Meta Estimate (Out State)' ~ 'Out of State'
  )) %>%
  ggplot(aes(x=shooting, y=Coef, ymin=lower, ymax=upper)) +
  geom_hline(yintercept=0) +
  geom_pointrange(fill='white', shape=21) +
  coord_flip() +  theme_bw() +
  theme(legend.position='bottom') +
  labs(x='', y='Meta Estimate Petitions')

#ggsave(petition_plot, width=8, height=5, 
#       file='figures/petition_state_RDITs.pdf')

ggsave(meta_plot, width=4, height=2, filename = 'datasets/petition_state_RDITs_meta.pdf')

